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Abstract 

With the question, "Can relativistic charged spheres form extremal black 
holes?" in mind, we investigate the properties of such spheres from a classical 
point of view. The investigation is carried out numerically by integrating the 
Oppenheimer-Volkov equation for relativistic charged fluid spheres and find- 
ing interior Reissner-Nordstrom solutions for these objects. We consider both 
constant density and adiabatic equations of state, as well as several possible 
charge distributions, and examine stability by both a normal mode and an en- 
ergy analysis. In all cases, the stability limit for these spheres lies between the 
extremal (Q = M) limit and the black hole limit (R = R+). That is, we find 
that charged spheres undergo gravitational collapse before they reach Q = M, 
suggesting that extremal Reissner-Nordtrom black holes produced by collapse 
are ruled out. A general proof of this statement would support a strong form 
of the cosmic censorship hypothesis, excluding not only stable naked singular- 
ities, but stable extremal black holes. The numerical results also indicate that 
although the interior mass-energy m(R) obeys the usual m/R < 4/9 stability 
limit for the Schwarzschild interior solution, the gravitational mass M does 
not. Indeed, the stability limit approaches R+ as Q — ► M. In the Appendix 
we also argue that Hawking radiation will not lead to an extremal Reissner- 
Nordstrom black hole. All our results are consistent with the third law of black 
hole dynamics, as currently understood. 
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1 Introduction 

Extremal black holes, black holes for which the charge or angular momentum param- 
eter equals the mass, occupy an exceptional position in black hole thermodynamics: 
due to their vanishing surface gravity they represent the absolute zero state of black 
hole physics. Over the past five or six years, compelling evidence has accumulated 
that extremal black holes represent a fundamentally different class of objects than 
their nonextremal counterparts. In some respects this is not surprising because the 
horizon structure changes completely at extremality, (see, eg. [|IJ]) and all the thermo- 
dynamic properties of black holes depend on the horizon structure. The exact nature 
of extremality is still the subject of some debate. Semi-classical calculations for eter- 
nal, extremal black holes || || indicate that the entropy is zero, while more recent 
semi-classical calculations for extremal objects collapsing into black holes ||, find 
that the temperature, and hence entropy, is not well defined. Similar results have 
been found for BTZ black holes |J. String theory calculations for extremal black 
holes, on the other hand , have yielded the ordinary Bekenstein-Hawking value for 
the entropy. 

Given the pathological properties of extremal black holes, at least from the semi- 
classical side, one reasonable supposition is that such objects are, for reasons as 
yet unclear, disallowed by nature. Israel |J, for example, has proven a form of the 
third law, which shows that it is impossible to reach extremality by any finite-time, 
continuous accretion process (with some caveats). On the other hand, some simple 
solutions of the Einstein equations for extremally charged objects are known || |10|) 
[TTJ . These studies, however, consider the rather unrealistic scenario of collapsing, 
infinitely thin shells. Moreover, the extremal solutions are evidently unstable, which 
again casts doubt on their relevance as models for genuine physical objects. 

In this paper, we take a purely classical point of view and consider the stability 
of charged, spherical matter distributions that satisfy the Einstein equations. Essen- 
tially we investigate interior Reissner-Nordstrom solutions, focusing on the extremal 
Q = M case. A study of relativistic charged spheres was undertaken thirty years 
ago by Bekenstein |T2|] , but without numerics he was unable to reach many firm con- 



clusions about the ability of charge to prohibit collapse. More recently a numerical 
investigation has been carried out by de Felice et al. JTB] (henceforth dFSY). The 
dFSY study was limited in that it considered only spheres with constant matter den- 
sity p and a power-law charge distribution. Furthermore, although p was taken to 
be constant in deriving the equation of hydrostatic equilibrium, the stability analysis 
assumed an adiabatic equation of state, i.e., p oc p^ m , where p rm is the rest mass 
density. As will become clear below, the p = constant case is the only convenient one 
for numerical integration; nevertheless, it is possible to construct a scheme capable of 
handling the more general adiabatic equation of state, as well as other charge distri- 
butions. We have done this for the present study and find that in all cases examined, 
instability sets in before extremality is reached. 
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2 Einstein and Oppenheimer-Volkov equations for 
charged spheres 

To derive the relativistic hydro- and electro-dynamic equations for a charged fluid 
sphere, we assume a general spherically symmetric metric in the form 

ds 2 = -e 2 *dt 2 + e 2A dr 2 + r 2 d6 2 + r 2 sin 2 £# 2 , (2.1) 

where $ = $o( r ) + t) and A = Ao(r) + Ai(r, t). In this section we are concerned 
only with the zeroeth-order (static) Einstein equations and so set $i = Ai = 0. The 
first-order quantities will be employed for the stability analysis in §^j. We write the 
Einstein equations in the form 

G^ v = R^ u — -g^R = 8nT^ u , (2.2) 

with c = G = 1. For charged spheres, the energy-momentum tensor will consist of a 
perfect-fluid part 

hydro 

= (p + p)u^u v + pg^ (2.3) 

plus an electromagnetic part 

MT^)em = F«F ua - -^F^F^. (2.4) 

In these expressions p = p rm + ^ is the total mass density, p rm is the rest mass density, 
e is the internal energy density, p is the fluid pressure, is the four-velocity, and 
F^ is the electromagnetic field strength tensor. As in the derivation of the exterior 
Reissner-Nordstrom solution, we take the magnetic field B to be zero. Due to spheri- 
cal symmetry, the electric field E must be entirely radial and the electromagnetic field 
tensor F^ u , which has only two components, F 01 = —E and F w = E, must satisfy 
the Maxwell equations 

1 -g-Fn,v=4nf, (2.5) 



where is the four-current. The only nonvanishing derivative for the static case is 
v = r. Working out ( |2.5p gives 

E(r) = £ - jJm, (2.6) 

where 

Q(r) = 4vr / e^ +A ^r' 2 f(r')dr' . (2.7) 
Jo 

Eq. (|2.6|) immediately implies 

F oi = e 2 W> . (2.8) 
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Note that ( |2.6| ) and ( |2.8| ) are of exactly the same form as the exterior Reissner- 
Nordstrom solution, as they must be by Gauss' law. The only difference is that for 
the interior case Q represents Q{r) rather than the total charge. {T^ u )em must also 
be of the standard Reissner-Nordstrom form, or from (|2.4j) and 



Q\r) 



(T^em = ^ dia 8[ e > - e > r > r sin 6 1- ( 2 - 9 ) 

Also note that T E m = (T^em — due to the antisymmetry of F^ v . 

The left-hand-side of the Einstein equations are the same as for any static, spher- 
ically symmetric solution. With the above expressions for T^, the (00) equation is 
found to be 

1 2A'e- 2A e- 2A Q 2 , N 

- + - = 8vrp+\, 2.10 

and the (11) equation is 



1 2$ ; e- 2A e' 2A _ Q 2 

ry 2 y^<2 

where (') denotes derivative with respect to r. 
If one multiplies ( |2.10|) by r 2 , one gets 



iTTp, (2.11) 



or, 

where 

m(r) 



-2.\ _ , -"'('•) F{r) 



r r 

Q 2 (r') 



An j p r' 2 dr' and .^""(r) = / ~ ) 2 ' dr' . 



In the case where p is a constant and Q is taken to obey a power law in r, these 
expressions reduce to those in dFSY. However, we will not restrict ourselves to this 
situation. 

To eliminate the metric functions $ and A and to get the basic hydrodynamic 
equations we consider the conservation laws T^ u v = 0. The fluid part is the same as 
can be found in the derivation of the ordinary Oppenheimer-Volkov equation (see, eg. 
14j|) and is found to be 

T- r = e- 2A b' + $'(p + P)]- (2.13) 
The electromagnetic part is found to be 

e~ 2A QQ' 



jrprr \ 
\ 1 EM);r 



4 vrr 4 
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Setting the sum of these terms to zero and solving for $' yields 



1 



P + P 



QQ' 



P 



Solving ([2.1 1|) for $' and inserting the result into (|2.14|) yields 



P 



P + P 



2r 



+ e 



2A 



Anrp + 



1_ Q! 

2r 2r 3 



47rr 4 (p + p)' 



(2.14) 



(2.15) 



With the form ( |2.12| ) for e 2 and some rearrangement of terms this can be written 

as 



P 



QQ' (p + p) 

4-7rr 4 r 2 



'„ 3 / A ^ 

Airr p + m(r) H 

^ v ; 2 2r 



2m(rl 



r 



(2.16) 



This is the Oppenheimer-Volkov (OV) equation for relativistic charged spheres, which 
will be the basis for our analysis. 

We point out that the Newtonian version of this equation is 



P 



QQ' pm(r) 
47rr 4 r 2 



which shows that in the Newtonian limit the Coulomb repulsion opposes the grav- 
itational force, as expected, helping to stabilize charged spheres against collapse. 
However in the relativistic OV equation with power-law charge distribution Q ~ r k , 
T/2 — Q 2 /(2r) < 0. Thus these terms decrease the effective mass in the numerator 
of the second term, tending to "stabilize." At the same time charge decreases the 
denominator, tending to "destabilize." The final outcome is not entirely clear, which 
is why one needs to investigate the OV equation numerically. 



Integration of the Oppenheimer-Volkov Equa- 
tion 



The basic strategy is simply to integrate ( |2.16|) and test for the stability of solutions. 



Before doing this one must specify the charge distribution and an equation of state. 
We have no serious models of any sort for charged "stars," or even know whether 
such objects exist. It seems intuitively reasonable that due to electrical repulsion the 
charge distribution should be weighted toward the surface, but otherwise all choices 
are basically ad hoc. One can imagine several plausible distributions. dFSY confined 
themselves to the simplest case, Q oc r k . We consider two distributions: 

Q(r) = Qt (^) fc e crtR = fPr k e cr / R , (3.1) 
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where R and (3 are adjustable input parameters that define the radius of the star and 
the total charge Qt respectively, and 



which flattens out for r 3> cR. Note that in both these cases the charge density, 
obtained by differentiating (|2.7| ), diverges at the origin unless k > 3. In point of fact, 
although we have been able to find solutions for the charge distribution ( |3.2|) , the 
stability analysis is generally inconclusive since no absolute boundary separation can 
be found consistently over all (3 of a given 7 across the entire range of parameters 
considered here. The one exception is 7 = 4 in which case the results are qualitatively 
similar to the exponential distribution ( |3.1|) for all (3, and so we confine ourselves to 
that case. 

Scarcely less arbitrariness applies to the equation of state (EOS) than to the 
charge distribution. Three obvious possibilities suggest themselves: p = constant, 
the adiabatic EOS p oc pl m with arbitrary adiabatic index 7, and the ultrarelativistic 
EOS, p = p/3. However, despite the fact that an exact solution with infinite central 
density p oc r~ 2 for the last case is known when Q = (see [0), we have found 
it difficult to find convergent solutions for this EOS when Q ^ when iterating 
over central pressure. (This is as opposed to iterating over central density, which 
we do for arbitrary adiabatic index, below). Thus we confine ourselves to the first 
two choices. We do, however, consider the 7 = 4/3 adiabatic EOS, which includes 
relativistic matter with non-negligible rest mass density (typically p r m/e ~ 10 2 — 10 4 
at the origin for 7 = 4/3, depending on p c and /3, and substantially lower ratios for 
larger 7). We reiterate that p represents the total mass-energy density of the matter 
plus gravitational field. Generally, the EOS is assumed to relate only the rest energy 
density, p rm , and the pressure. As mentioned in §|2|, the total local energy density is 
P = Prm + e, where e is the internal energy density. We consider only polytropes, for 
which e = p/(j — 1). 

Once the charge distribution and EOS are specified we can integrate the OV equa- 
tion. For virtually any case one can imagine, however, integrating ( |2.16| ) is impossible 
analytically and the integration must be carried out numerically. Even then, except 
for the case p = constant, this turns out to be nontrivial. For p = constant, the 
total mass will be simply m(r) = (4/3)irr 3 p. With the boundary condition p(R) = 
and a form of Q to make T analytically integrable, this allows convenient integra- 
tion from the outer boundary inward. For more reasonable equations of state, how- 
ever, one expects p to vanish at r = R as well (eg., if p = kp 1 ^" 1 +p/('j — 1)). In 
such cases the entire second term in ( |2.16| ) vanishes at the outer boundary, which 
causes a serious problem numerically because dp/dr\ r= R = QQ' j (Airr 4 ) > 0, while 
p(R — Ar) ps — QQ'\ Ar|/(47ri? 4 ). Thus, unless QQ' < 0, an inward integration 
scheme immediately predicts a negative pressure unless the integration is initiated 
from some arbitrary outer pressure and surface radius values to compensate for the 
charge distribution. 




(3.2) 
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For this reason, we choose to integrate outwards from r = 0. This strategy also 
presents difficulties. We can, say, choose R as the prespecified "radius" of the star, 
but in starting the integration at the origin, we have no prior knowledge of which 
combination of parameters p c (central density), p c (central pressure), 7, /3, k, c, 
and Q will actually produce a star with radius R. In other words, if we define the 
physical surface radius R p to be the point at which the pressure first goes negative, 
there is no reason that R p will equal R. We have therefore adopted the following 
procedure: For a given parameter set (R, /3,j,k,c,Q, p c ,p c ) we take a first guess at 
p c = 1/^Rp]/ 1 ) with arbitrary p c , then iterate on p c using a bisection algorithm until 
R p converges to R. The bisection method assumes a locally monotonic correlation 
between 5p c /p c and 5R/R, where the (n + l)-th estimate of the central density is 



given by p c n+1) = p c n) + 5p c n) . One chooses 5p{ n) by p c n) 5R/R = ap c n) (R p - R)/R, 



where a is a coefficient initially set to unity, but is halved for convergence whenever 
SR changes sign. We find in most cases we can converge on R to a relative tolerance 



The numerical integration is performed using an adaptive fourth order Runge- 
Kutta scheme to maintain a constant error tolerance of 10 -12 over a grid constructed 
with cell spacing scaled to the surface radius such that Ar/R < 5 x 1CT 4 . Further 
accuracy is achieved at the outer surface where pressure and density vary sharply 
by adaptively refining the mesh spacing to a minimum level Ar « 10~ 12 when the 
pressure is observed to change sign. This procedure works well if R is sufficiently 
small. However, as R is increased one eventually hits a barrier R = R max beyond 
which convergence to the 0.1% tolerance is impossible or (2m(r) + F)jr becomes 
so large that the metric function (|2.12|) becomes negative. This latter marks the 
black-hole limit R = R + = M + a/M 2 — Q 2 , where M is the exterior gravitational 
mass. The upper bound on R max generally depends on all input parameters, and is 
computed numerically for each case presented in this paper. For example, considering 
an adiabatic EOS with Q = \^j3r 3 , Figure p] shows the maximum surface radius as 
a function of (3 for the different 7 considered in this paper. We are unable to find 
convergent solutions at radii above each of the curves. 

Figure || plots the total mass m(R) as a function of surface radius for 7 = 4, 100 
and fixed (3 = 1. Here the surface radius is increased from 5% of R m ax to the maximum 
computed value R ma x- Each of the solid lines represent different central pressures p c 
with the larger pressures resulting in greater masses. Also plotted in Figure |^ are the 
corresponding black hole and extremal charge solutions for comparison. The black 
hole solution for r = R is found by setting the metric function ( |2.12| ) equal to zero at 
that point, which gives m(R) = (R—J r (R))/2. By contrast, the extremal solution (not 
necessarily a black hole) is found by first matching the exterior Reissner-Nordstrom 
metric to the interior solution (|2.1|) at r = R: 



of 5p/p < 0.1%. 




(3.3) 
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which gives the definition of gravitational mass 

1 f R , 9 Q 



1 r R O 2 O 2 

M = U i^pr 2 + ^)dr + ^ 5 . (3.4) 



Setting Q = M gives the extremal curve for m(R) = Arc J pr dr: 

m{R) =Q~^-^T- (3-5) 

Finally we note that the total charge to total mass ratio Q/M generally increases 
monotonically with increasing R and fixed /3, as Figure [3] demonstrates. Since we are 
interested in the maximum or extremal limit of Q/M, the subsequent stability plots 
presented in section f| are computed for those parameter values at the upper bounds 
of R (= R m ax) for which we are able to find convergent solutions. 

As a reminder, we point out that because we do not know beforehand the exterior 
mass of the star, the radius R in our graphs is not given in units of M. R is merely 
where the pressure drops to zero expressed in geometric units. 



4 Stability I: Normal Mode Analysis 

We provide two stability analyses. The first, discussed in detail in this section, is 
a normal-mode analysis for radial pulsations of the charged sphere, similar to that 
performed by dFYS (see also or ]TB[). The second, discussed in §EL is based on 
the variation of energy. 

We consider radial pulsations of the sphere with the goal of writing the Sturm- 
Liouville equation: 

(FC)' + (H + cu 2 W)C = 0. (4.1) 

Here ( is a "renormalized displacement function," and F, H and W are functions of 
the zeroeth-order (equilibrium) solution to the Einstein equations, to is the angular 
frequency of the assumed sinusoidal time dependence of Q. Once these functions are 
known, we evaluate 

and deduce that the solution is stable under radial perturbations if uo 2 > for ap- 
proximate trial functions £ that satisfy the same boundary conditions as (. 

To derive F, H, and W is a tedious undertaking; because we are now dealing 
with a time- dependent problem, the metric functions in (|2.1| ) must be taken to be 
$ = $o( r ) + &i{ r ,t) and A = A (r) + Ai(r, t), where we assume the first-order 
quantities are small. Once the perturbation equations are derived, to get the Sturm- 
Liouville form, the first-order quantities are eliminated in favor of the zeroeth-order 
solutions. Nevertheless, because our equations differ slightly from dFYS, we now 
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outline the procedure. The development basically follows that of MTW JTSJ, chapter 
26. 

We first consider the total energy-momentum tensor. With the condition g flu u fJ 'u u = 
— 1, and raising and lowering with the full metric, one gets to first order w = 
— e*°(l + $i) and U\ = e 2A ° - * ^, where £(r,t) <C r represents the small displace- 
ment of a fluid element initially at radius r. For the electromagnetic part of T^ u , 
we let Q = Qo + Qi, with Qi = —Q' (r)£(r,t) to first order. This is equivalent to 
assuming that the charge distribution in the oscillating system is the same as in the 
unperturbed system, i.e, that the Lagrangian perturbations vanish and no electric 
currents are introduced for the comoving observer. Then, 

2*0 i 2#n I i / Qo i o "\ QoQ'o£\ 

T °° " [ Po + 8^) e +e l v Pl + $l( 4^ + 2po) ~l^J 
T i = -e 2Ao (Po + Po)e 

Hi = L- A^ 2Ao +e 2Ao L+A 1 (2 Po --^-) + %^V (4.3) 



4 / \*"- ^ 47rr 4; 47rr 4 



Note the presence of T i, due to the fluid motion. This is entirely a first-order quantity. 

The left-hand-side of the Einstein equations are most easily computed with the 
help of Mathematica or Maple. Then the (00) Einstein equation is found to be 



e 2 * , 2 



1 + 2$!) + -e" 



-2A +2<I>o 



^"2 iy 2 



- + rA; + Ai(l - 2rA;j -*!(!- 2rtt ) + rA[ 



Equating orders on each side leads to the O t/l -order equation ( |2.10| ), which minus 
the charge term is also MTW Eq. (26.1b). The first-order equation, after some 
simplification, comes out to be 

[A x (l - 2rA' ) + rA'J = 8vr Pl - ?9^.. (4.5) 
With Tqi given by ( |4.3| ), the (01) Einstein equation is then nontrivial and gives 
Ai = -47rre 2A "(po + Po)£ = ~(K + *d)£, ( 4 -6) 



where the second equality follows from the zeroeth-order equations ( |2.10 ) and ( 2.11|) . 



We will shortly need an expression for the perturbed pressure p\. One can get this 
from the law of baryon conservation d(An)/dr = —n(V • u) = (—g)~ l l 2 (^/—gu a ), a . 
Taking a — t,r yields 



An = -w r-V A °(r 2 e Ao O' + Ai , (4.7) 
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where An represents the Lagrangian perturbation in n. Assuming adiabatic fluctua- 
tions such that Ap/An = jp/n then gives 



Pi 



-IPo 



r -2 e -Ao/V2 c Ao£V 



r 2 e Ao O' + Ai -&>o 



which is MTW Eq. (26.9). 

The first-order (11) Einstein equation is found to be 



=4vrre 2Ao [p 1 + 2A 1 ( y p 



Qo 



, A i e2A ° , e 2A gogoe 



r 



(4.9) 



With ( |4.6|) and (|4.8|), Ai and p\ can be eliminated to get 



-4vrre 2A ° 
_ 47rre 2A ^ 



7P r~V At, (r 2 e Ao £)' 
j/ + (A' + $' ) ( 2 Po - JPo -M 7+ ' 



47rr 4 47rr 2 



4 vrr 4 



(4.10) 
, = 0, 



By considering the parallel component of the conservation equation u^T^ 
using ( P~7| ) and assuming the Lagrangian perturbation AQ vanishes, one can show 



Pi 



-(Po + Po) 



-V A °(r 2 e A °£)' + A a -e^b, 



(4.11) 



which is MTW Eq. (26.11). 

The transverse component of the conservation equation, h^T a ^. a = 0, where the 
projection tensor h% = + u^u v eventually yields 



3 2A -2$o 



(Po + Po)£ 



%(Pi+Pi) -^[(Po+Po) -p'l 

^(Qfe + QoQot + QoQ'o? 



4nr 4 



(4.12) 



We note that the Qq term does not appear in dFSY. 

With these expressions we are finally able to write down the Sturm-Liouville 
equation. We define the renormalized displacement function ( = r 2 e~*°£ and assume 
( = ((r)e~ iujt . Then 



F 

W 

H 



7Po^" 2 e Ao+3$0 , 



r- 2 (p +Po)e 3Ao+,I>0 , 

4vrr(p + p )e 2A ° (-p' + ^ + ±(p + Po] 



r -2 e A +3$ 



(4.13) 



+2p^' + pg + $^ _ M _ ^(q* + g gg + g g(,($' _ 2)) 



We have not attempted to reconcile this expression for H with the one in dFSY. 

Physically reasonable solutions require (/r 3 finite or zero as r — > and 7Po r_2 e*°C / ' 
as r — > i?. Any trial function £ that satisfies these conditions and that extremizes 
(4.2) is acceptable. Following Chandrasekhar and dFSY, we choose £ oc r 3 . 
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Figure ^ shows the stability curves for the adiabatic EOS with Q = y^/Sr 3 . The 
results are generated by choosing an initial central pressure (equal to unity) and 
iteratively solving the OV equation to convergence in the central density for fixed j3 
and surface radius R = R ma x- The central pressure is then monotonically increased 
and the OV equation repeatedly solved until the angular frequency defined in equation 
( |4.2| ) becomes negative. The data points connected by solid lines are the solutions 
at which the frequency first becomes negative. Also shown are the black hole limit 
R = R + and the Q = M limit, generated from Eq. fl3.5|) . Note that the stability limit 
falls between the Q = M curve and the black-hole curve. As can be seen from Figure 
[|, Q/M varies over a broad range at R = R ma x, but is constrained by Q/M < 1 in 
the limit of large 7. Thus instability sets in before extremality is reached since Q/M 
decreases as the R = R + curve is approached (say by increasing the central pressure) , 
even in the limit 7 — > 00 where all the curves come together. 

Although for reasons of clarity we have not shown the 7 = 4/3 and 5/3 solutions 
in Figure |], they are qualitatively similar to the other cases, except that the extremal, 
stability and black hole curves are more widely separated. A more detailed plot of the 
7 = 100 case is shown in Figure [| for two different values of j3, and varying the surface 
radius up to R = R ma x- Also shown in Figure |5] with solid lines are the corresponding 
limits for the interior Schwarzschild solution when Q — 0, m(R) = 4R/9. This 
absolute limit is a good match to the numerically computed results even for relatively 
large ratios of Q/M and radii R ~ R ma x- However, a more useful extension of this 
limit should employ the gravitational mass M deriving this constraint. Figure |6] plots 
the same result as|5|, but for the gravitational mass M(R). The solid lines are now the 
limits for the exterior Schwarzschild solution M(R) = 4R/9. Notice that the upper 
bound of stability with this definition is always below, but approaches the black hole 
limit R — > R + — > M as Q — > M and R —>■ R ma x- This result is consistent with that 
of [|l7j who find no non-singular static solutions with R —>■ R + for Q < M. 

The overall behavior is repeated in Figures ^ and || which show similar graphs for 
two additional cases: the adiabatic EOS with Q = y^j3r 3 e r ^ R , and a constant density 
profile respectively. Once again, unstable solutions separate the Q = M and black 
hole states. This is our main result: that the onset of instability always occurs before 
extremality is reached. 



5 Stability II: Energy Analysis 

Admittedly, the above stability analysis is not exceptionally transparent. Somewhat 
more insight into the stability of relativistic charged spheres can be gained by consid- 



ering an energy analysis. Following [16] and neglecting the electrostatic energy which 



is invariant with respect to changes in density, one can write the internal (thermal 
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plus gravitational) energy of a spherical gas cloud as 



E 



P 1 



2m(r) 



r 



1/2 



Pr 



dV, 



(5.1) 



where the invariant volume element dV = 47r(l — 2m(r)/r — J r (r)/r)^ 1 ^ 2 r 2 dr. An 
equilibrium configuration is determined by the condition dE/dp rmc = 0, where p rmc 
is the central rest-mass density, and the transition from a stable to unstable configura- 
tion is determined by the condition d 2 E/dp 2 mc = at the point where dE/dp rmc = 0. 
We find the transition numerically by first iteratively solving the OV equation for a 
given central pressure p c until convergence is reached in p c . The energy integral (|54 



is then evaluated and tabulated as a function of total mass and central rest-mass den- 
sity. The inflection point is then easily computed as a function of M by discretizing 
the tabulated data with respect to p rmc . Figure |] shows a comparison of this method 
of evaluating stability with the normal mode analysis of §^ for the 7 = 4 case of 
Figure [|, but on a linear-linear scale. The results agree remakably well. However, 
we note that this energy approach is much more sensitive to numerical resolution 
and accuracy than the normal mode method. In fact, this method proves generally 
inconclusive for extreme (both small and large) values of 7. Hence we rely exclusively 
to the more robust normal mode method of §(| to assess stability. 

Eq. ( |5.1|) is already suggestive in that one can immediately see that the function 
T enters with the same sign as mir). In other words, charge effectively increases the 
magnitude of the gravitational potential energy, against the internal energy included 
in p, which should tend to destabilize the charged sphere. This can be confirmed by 
a perturbation analysis. Let p = p rm (l + u), where now u is the internal energy per 
unit mass, and assume that both k C 1 and m(r)/r <C 1. Note that J- '/r ~ Q 2 /r 2 = 
(Q 2 1 m 2 ) (m 2 1 r 2 ) , so this is at least a second order effect, T jr ~ m 2 /r 2 if Q ~ m, and 
even higher order if Q <ti m. Expanding (|5.1|) to lowest order in J 7 , gives 



E 



R 



Pr 



m 



u — 



urn 
r 



m 



2r 



dV. 



(5.2) 



The first four terms are those found in computing the standard relativistic cor- 
rections to a Newtonian star (see flrEfl , chap. 6). The last term, which we denote as 
Iq = (1/2) / p Tm [T jr)dV is the correction due to charge. Note that Iq > always. 
For a 7 = 4/3 polytrope one finds with the correction due to charge 



E — {AM — BM^)p\t + CMp-^ - DM 1 " p 2 r t - I Q , 



(5.3) 



where A, B, C, D are positive numbers that depend on fundamental constants as well 
as on results of the numerical integration (see |16|, Eq. (6.10.23)). 
The condition dE/dp rmc = gives 



\{AM - BM*l*)p££ - \CM P -^ 



DM 7 " n^ 1 " - V 



0. 



(5.4) 
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Here, I'q = OIq/ 'dp rmc ~ (1/2) /(a(r)jF '/r)dV ', where a(r) represents the normalized 
density profile p rm = p rmc a(r) with a(r = 0) = 1. We have also neglected higher order 
contributions from the invariant volume element dV. For real stars, equating the first 
term in equation ( |5.4|) to zero leads to the Chandrasekhar limit. The remaining terms 
are small corrections. Solving for the first term in the above equation and using the 
result in the condition d 2 E/dp 2 rnc = yields 

\cMp~U 3 - 2 -DM 7 ^p-^ - \p~Llq = 0. (5.5) 

To leading order M ~ (A/B) 3 ^ 2 , so the onset of instability occurs at 

C(^ /2 - D{^yi 2 Pcrit - 3p^ t I' Q = 0. (5.6) 

When Iq = 0, this gives the standard result p crit = (C/D)(B/A) 2 . However, since 
I'q > 0, the presence of the last term clearly lowers p cr u, decreasing the threshold of 
instability by making the star unstable at lower densities and masses m(r) at fixed R. 
This behavior is observed in Figure |5] where as Q and Q/M are increased the stability 
curve falls further below the R = (9/4)m(R) Schwarzschild stability limit. But, since 
we are increasing charge, we are increasing the gravitational mass M over m(r) (see 



O]) , and so, paradoxically, when measured in M, the stability limit increases, as seen 



in Figure 



6 Conclusions 

We have seen that, contrary to Newtonian intuition, in general relativity charge does 
not tend to stabilize or counterbalance the interior energy density of fluid spheres. 
To be sure, our numerical results indicate that for a fixed radius, instability sets in 
at a smaller mass m(r) than in the fluid only case. Nevertheless, charge does in- 
crease the gravitational mass M, and the instability threshold for relativistic charged 
spheres monotonically approaches the black hole horizon in the extremal limit. Most 
important, though, for all cases studied the onset of instability of relativistic charged 
spheres takes place before the attainment of extremality. At that point, gravitational 
collapse into black holes takes place, but these are not extremal holes. All these 
results are certainly consistent with the currently accepted form of the third law of 
black hole dynamics, which forbids attainment of the extremal state by finite time 
processes involving accretion of positive energy f| (see also Appendix). 

One might legitimately doubt whether studies such as this one have any bear- 
ing on reality given that astrophysical objects tend not to be charged. For precisely 
this reason it has long been doubted that Reissner-Nordstrom solutions can repre- 
sent astrophysical objects. On the other hand, to the extent that charged solutions 
shed light on the third law, they have been of great interest. Regardless, however, of 
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whether one chooses to interpret charged spheres as stars or as solitons or as mod- 
els for elementary particles, our results make it difficult to see how extremal objects 
could arise from any remotely realistic collapse scenario. A general theorem ruling 
out stable extremal black holes altogether would be desirable; it would effectively be 
an alternate formulation of the third law, as well as a strong version of the cosmic 
censorship hypothesis (or theorem, in that case); not only are stable naked singular- 
ities excluded, but stable extremal black holes as well. 
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Appendix: A note on Hawking radiation and the 
third law 

Israel's proof of the third law || , that the surface gravity of a black hole can never be 
forced to zero in a finite-time, continuous process, assumed accreting charged shells 
and that the energy flux crossing the outer apparent horizon of the black hole was 
positive. Therefore the proof has nothing to say about Hawking radiation, because in 
this case a negative energy flux crosses the apparent horizon, violating the assump- 
tions of the theorem. Moreover, at first glance, it might seem that Hawking radiation 
would drive a black hole toward extremality: If a black hole is emitting neutral parti- 
cles, for example, the charge would remain constant while the mass decreased, forcing 
Q/M — > 1. It is not terribly difficult to see that this does not happen. 
The power radiated by the black hole is 

P = ^~R*T A ~R i K*>, (A.I) 
dt 

where k is the surface gravity. For a Schwarzschild black hole, R = 2M and k = 
1/4M. This immediately gives the well known result for the evaporation time r ~ M 3 . 

For a Reissner-Nordstrom black hole, however, R = R + = M + y/M' 2 — Q 2 and 
the surface gravity is given by 

(M*- Q r 2 (A2) 

2M[M + (M 2 — Q 2 ) 1 / 2 } — Q 2 y ' ; 

(see Wald |L8| Eq. (12.5.4)). Of course, as Q — > M, k approaches zero, and therefore 
the flux, which goes as the fourth power of k is rapidly decreasing. This suggests that 
an infinite amount of time is required to reach Q = M. The statement can be made 
rigorous by integrating ( |A.1|) for Reissner-Nordstrom black holes. With R = R + and 
the above expression for k we have 

f (2M(M + VM 2 -Q 2 )-Qr 

J (M + ^M 2 - Q 2 ) 2 (M 2 - Q 2 ) 2 ' 1 j 
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where the limits should be from M to Q, assuming M > Q initially. With Q = 
constant, the expression can be integrated analytically. The result for the indefinite 
integral is 

32M 3 ( ie 2 , MQ 4 



r (M) = — — + 16MQ 2 + 



3 ^ 2(Q 2 -M 2 ) 

2(16M 4 + 16M 2 g 2 - 41Q 4 ) 35Q 3 tanlr^M/Q) 
+ 3^M 2 - Q 2 2 ' ( } 

Assuming, for example, that M initia i ^> Q and Mf ina i = Q(l + e) with e < 1, then 

64M 3 g 3 ^. _ 1/2 . . 

T — + ^ + 0(e 1 ' 2 ), (A.5) 

which clearly diverges in the limit Q —>■ M (equivalently e — > 0). 

Although we have done this only for Q = constant, it may hold true generally. 
If the Hawking temperature T is initially below the rest mass of the electron, the 
black hole will radiate only neutral particles and so Q will remain constant while 
M decreases. Plotting ( |A.2| ) shows that k and hence T will actually increase with 
decreasing mass until Q = (a/3/2)M, at which point k = 2/(9M) then monotonically 
decrease to zero as n ~ \/2e/Q in the limit M = Q(l + e). Assuming the maximum 
temperature is always less than the mass of the electron, then Q can be expected to 
remain constant, and evaporation cannot proceed in finite time, as discussed above. 

If at any time T becomes greater than m e , the hole will evaporate charge. But 
for any charged elementary particle g>m due to negligible size of the gravitational 
interaction. Assuming that for temperatures above the rest mass, all species of par- 
ticles are radiated with roughly equal probability (equipartition), under evaporation 
Q will become negligible compared to M and extremality will never be approached. 

We emphasize that, as Israel points out, the entire discussion neglects any super- 
radiant component which is nonthermal, but it is not clear whether this would alter 
the conclusion. We have also neglected other processes, such as discharge from pair 
production mentioned by Page JTP|, which in this case would likely only strengthen 



the conclusion. In any case, certainly for Q = constant, Hawking radiation cannot 
violate the third law of black hole dynamics and perhaps not in general. This is 
consistent with the early work of Page and with the point of view that extremal black 
holes do not exist. 
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Figure 1: Maximum surface radius as a function of charge parameter f3 found 
by converging in the central density to a tolerance of Ap c / p c = 1CT 3 for a charge 
distribution Q = and the different adiabatic indexes considered in this paper. 

For fixed /3, the numerical solutions do not converge to the specified tolerance at radii 
above the plotted curves. 
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Figure 2: Total mass density m as a function of outer surface radius i? (up to 
R = Rmax) for fixed charge parameter (3=1. Each of the solid lines are solutions 
to the OV equations with different central pressures (the larger pressures correspond 
to greater mass densities). Also plotted are the densities for the corresponding black 
hole (R = R + ) and extremal charge (Q = M) limits. 
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Figure 3: Total charge to mass ratio Q/M is a monotonic function of surface radius 
up to the maximum computed radius R max , and peaks at Q ~ M for large 7. The 
results are shown for a central pressure p c — 1. 
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Figure 4: Total mass function corresponding to the stability limit as a function of 
the maximum computed surface radius for 7 = 4, 10, 100. The solid lines are found 
by increasing the central pressure from unity and plotting the first points which 
become unstable. The four points in the 7 = 4 curve represent for this distribution 
Q = v^r 3 , = 0.01, 0.1, 1.0, 10 (right to left). The three points for 7 = 10 represent 
P = 0.01, 0.1, 1.0; the two points for 7 = 100 represent (3 = 0.01, 0.1. The displayed 
data are chosen to prevent overlap of the different 7 curves, but we note that the 
general conclusion remains the same for all cases and choices of f3 investigated: the 
limit of stability lies between the Q = M curves (dashed lines) and the black hole 
R = R + limit (dotted lines), even for 7 = 4/3 and 5/3 which are not shown here for 
clarity. Thus instability sets in before extremality is reached. 



Instability of . 



21 



0.8 




0.2 1 1 1 1 1 1 1 

0.5 1 1.5 2 

R (y=100) 

Figure 5: As Figure |], except here we show closeups of the extremal, black hole, and 
stability curves for the more interesting high adiabatic index (7 = 100) case over a 
range of surface radii leading up to the maximum R = R ma x- Also plotted is the solid 
line m(R) = 4R/9. The usual stability limit for the interior Schwarzschild solution is 
actually given by M = 4R/9, where M is the gravitational mass. The two expressions 
are equivalent only in the limit Q = 0; see also Figure |6|. In all cases, the stability 
curves separate the extremal charge and black hole states. 
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Figure 6: As Figure [5], except here we plot solutions for the gravitational mass M 
as a function of surface radius R. This time the solid lines show the stability limit 
M(R) = 4R/9 for the interior Schwarzschild solution. Notice there is no absolute 



upper bound on stability other than the black hole limit R 
Q^M. 
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Figure 7: As Figure |], but for a charge distribution of Q = y/]3r 3 e r ^ R . The 
represent (3 = 0.01 (right-most points) and 1.0 (left-most points). 
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Figure 8: As Figure f|, but for a constant density distribution. The symbols represent 
(3 = 0.01, 0.1, 1.0, 10 (right to left). Although difficult to see, an unstable solution 
(sold line) always lies between the extremal and black hole curves for all surface radii 
R < Rmax we have studied. This figure is shown for R = R ma x as are most of the 
other stability plots. 



Instability of . 



25 



0.07 



0.06 






- y=4, normal mode 






G 


-o y=4, energy integral 


0.05 






- Q=M 








R=R 

+ 



0.04 - 



§• 0.03 - 
0.02 - 

0.01 - /^x^ 
0.00 

-0.01 ' 1 1 1 1 1 1 1 1 1 1 1 1 1 — 

0.00 0.02 0.04 0.06 0.08 0.10 0.12 0.14 

R 

Figure 9: Comparison of the energy variation and normal mode methods in evaluating 
the stability of solutions for the same 7 = 4 case as Figure |, but displayed on a 
linear-linear scale to distinguish the stability lines. Over the parameter range that 
both methods can be applied to reliably (mostly restricted by the energy method 
which is less robust from a numerical point of view), the results agree remarkably 
well. 



